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The  Chemical  Equilibrium  Problem: 
an  Application  of  SUMT 


MCLEAN,  VIRGINIA 


FOREWORD 


This  paper  deals  with  the  determination  of  the  equilibrium 
concentration  of  multiphase- reacting  chemical  systems.  This 
type  of  problem  is  of  paramount  interest  in  many  contemporary 
technologies  and  in  recent  years  has  attracted  the  attention  of  the 
mathematical  community.  The  reason  this  particular  problem  is 
of  interest  is  that  it  can  be  formulated  as  a  mathematical  non¬ 
linear  programming  problem.  In  keeping  with  RAC’s  research 
interest  in  nonlinear  programming,  in  particular  the  sequential 
unconstrained  minimization  technique  (SUMT),  it  is  felt  that  it  is 
of  great  importance  to  illustrate  the  power  of  SUMT  in  obtaining 
the  solution  of  moderate- sized  problems.  Inasmuch  as  many 
groups  are  concerned  with  developing  special  algorithms  for  solv¬ 
ing  the  chemical  equilibrium  problem,  the  application  of  the  gen¬ 
eral  SUMT  to  this  problem  is  presented  here. 

It  is  felt  that  it  is  necessary  to  document  the  practical  appli¬ 
cation  of  SUMT  as  well  as  the  theoretical  developments,  and  this 
problem  is  of  a  size  that  again  demonstrates  the  power  of  the 
technique. 


Nicholas  M.  Smith 

Head,  Advanced  Research  Department 
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The  Chemical  Equilibrium  Problem: 
an  Application  of  SUMT 


ABSTRACT 


The  equilibrium  composition  of  homogeneous  chemical 
reactions  is  a  very  important  problem  in  studies  of  complex 
chemical  systems  that  arise  in  the  study  of  certain  physio¬ 
logical  systems,  rocket-propulsion  systems,  reentry  prob¬ 
lems,  etc.  The  aim  of  this  report  is  to  show  that  the  sequen¬ 
tial  unconstrained  minimization  technique  (SUMT)  as  exploited 
by  Fiacco  and  McCormick  is  capable  of  solving  this  problem 
and  indeed  offers  distinct  advantages  over  the  heretofore  best 
available  computational  techniques  such  as  those  presented 
by  White  et  al1  and  DeHaven  and  Deland.2  The  power  of  this 
method  for  solving  moderate-sized  problems  is  illustrated 
by  the  choice  of  a  45-variable  problem  with  16  linear  equality 
constraints. 


1.  INTRODUCTION 

The  problem  of  determining  the  equilibrium  composition  of  homogeneous 
chemical  reactions  can  be  shown,  in  a  great  number  of  cases,  to  be  equivalent 
to  finding  that  composition  vector  that  results  in  minimizing  the  Gibbs- Helm¬ 
holtz  free  energy  of  a  chemical  system  and  simultaneously  satisfying  the  re¬ 
quirements  of  mass  balance.  This  observation  replaces  a  chemical  problem 
with  a  mathematical  problem:  find  a  vector  that  minimizes  a  certain  function 
subject  to  the  satisfaction  of  a  system  of  linear  equations  that  are  the  formal 
expression  of  the  mass-balance  law. 

The  .atent  of  this  report  is  to  show  that  SUMT  as  exploited  by  Fiacco  and 
McCormick3'6  is  a  method  that  offers  distinct  advantages  over  preexisting 
methods  that  were  developed  for  this  particular  problem.1’7  Until  recently,8 
all  these  methods  depended  for  their  solution  on  the  availability  of  positive 
estimates  of  the  solution.  Advantage  has  also  been  taken  of  the  special  structure 
of  this  problem  to  greatly  enhance  the  rate  of  convergence  of  the  minimization 
technique  utilized  in  SUMT. 

2.  THE  SUMT  ALGORITHM 

A  brief  discussion  of  SUMT  for  nonlinear  programming  is  presentee  in 
this  section.  No  pretense  is  made  of  being  very  detailed,  and  the  reader  is 
referred  to  the  basic  works  on  the  subject  as  given  in  Fiacco  and  McCormick.3-6 
Consider  now  the  solution  of  the  general  nonlinear  programming  problem: 

Minimize 

fix) 

subject  to  (A) 

h,  (x)  =  0,  j  =  1 . p 

8j  M  >0,  i  =  1 . m, 
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where  x  6  E",  i.e., 


The  SUMT  algorithm  for  the  solution  of  this  problem  is  as  follows: 

(a)  Form  the  function 

P(x,  rk)  ^  f(x)  +  r„  1  l/8j  (i)  +  rj*  I  hf  (x). 

(b)  Find  the  unconstrained  minimum  of  the  P  function  in  the  region 

t1  Ifli  (J)  >  0,  i  =  1,  .  .  . ,  m}  for  each  fixed  r^,  where  [ rj, )  is  a  decreasing  se¬ 
quence  of  positive  numbers  such  that  lim  rj,  =  0. 

k 

It  will  follow  that,  under  certain  reasonable  restrictions,  the  sequence  of 
values  of  the  P  function  {P(x  ;rk  )],  respectively  minimized  by  Cx  (f),) }  over  the 
strictly  monotone  sequence  { rk  ] ,  converges  to  the  optimum  value  of  f. 

At  this  point  it  would  not  be  amiss  to  point  out  that  the  method  SUMT 
uses  to  find  the  unconstrained  minimum  of  P(x;r)  is  a  modified  Newton’s 
method.  That  is  to  say,  starting  with  an  initial  estimate,  for  fixed  r,  one  de¬ 
termines  the  next  move  as  follows: 


x"  +  l  =  x"  -  0l72P(x";r;]-1  VP(*V). 


The  parameter  6  is  a  scale  parameter  whose  determination  may  vary  lending 
rise  to  different  types  of  algorithms  for  minimization.  This  point  is  not  dis¬ 
cussed  in  this  paper.  The  essential  thing  to  observe  here  is  the  difficulty  that 
could  arise  in  the  determination  of  the  inverse  of  the  matrix  v2P(x;r).  Note 
also  that  if  the  function  f  is  convex  and  the  functions  a,  are  concave  and  if  the 
hj’s  are  linear  then  the  resultant  P  function  is  convex,  and  hence  when  the 
feasible  region  has  a  nonempty  interior  the  problem  before  us  is  to  find  the 
unconstrained  minimum  of  a  convex  function. 

If  it  should  turn  out  that  the  matrix  v2P(x;r)  can  be  put  in  the  form  A  + 
uav1,  where  A  is  an  n  x  n  matrix,  u  is  an  n  x  p  matrix,  cr  is  a  p  x  p  matrix, 
and  vT  is  a  p  x  n  matrix,  then  the  following  formula  may  be  used  for  computing 
[  A  +  utn>T]-1. 

IA  +  uwT]  1  =  A~’  -  A-1u[<7~’  +  uT  A-1u]-1  uT  A-1.  (1) 
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r-V 


!■' 


If  the  matrices  A  and  cr  are  easily  invertible  then  it  is  obvious  that  this  method 
(“rank  annihilation”)  is  extremely  powerful. 

It  is  also  the  intent  of  this  paper  to  show  that  by  using  this  method  of 
matrix  inversion  the  SUMT  can  be  brought  to  bear  to  solve  some  medium- sized 
mathematical  programming  problems,  collectively  called  the  “chemical  equi¬ 
librium  problem,”  that  arise  in  various  scientific  and  technological  areas.  The 
main  advantages  of  this  method  of  matrix  inversion  when  applicable  is  that  it 
greatly  enhances  the  rate  of  convergence  of  the  SUMT  algorithm. 


3.  THE  CHEMICAL  EQUILIBRIUM  PROBLEM 


The  chemical  equilibrium  problem  is  the  determination  of  the  composition 
of  the  chemical  species  that  minimizes  the  (Gibbs)  free  energy  of  a  chemical 
system  and  must  also  satisfy  the  mass  conservation  law  in  the  form  of  the 
mass  balance  equations.  Mathematically  speaking,  the  problem  is  as  follows: 

Minimize 

F(x), 

where 

?(x)  =  kti  {,!,  Ifk  [Cfl>  +  f"  (v  *1  ^k)]} 


subject  to 
where 


Hr  =  b  and  r  >  0 


b  =  (bl . b,)T-  1  -  (*11.*2I . V.p)1 


In  this  description  a  chemical  system  (mixture)  is  considered  to  be  a 
system  composed  of  a  number  of  different  molecules  of  various  species  exist¬ 
ing  in  p  various  phases,  e.g.,  physical  compartments,  conceptual  compartments, 
or  a  gas,  liquid,  or  solid.  It  is  to  be  emphasized  that  the  same  molecule  in  a 
different  phase  is  to  be  considered  a  different  molecular  species.  Assume  that 
all  the  different  species  occurring  in  the  system  are  composed  of  some  finite 
set  of  basic  ingredients,  which  may  be  molecules,  atoms  (both  charged  and 
uncharged),  aud  possibly  pure  charge.  The  only  requirement  is  that  all  the 
reaction  products  be  combinations  of  these  basic  ingredients,  or  this  basic 
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“alphabet.”  The  determination  of  what  is  a  possible  reaction  product  is  of 
course  left  in  the  hands  of  the  chemist  who  might  be  applying  this  method. 
Presumably  he  can  determine  a  priori  what  the  possible  nontrace  amounts  of 
the  reaction  products  will  be.  For  instance,  ir  there  are  bj  atomic  weights  of 
basic  unit  j  originally  in  the  system  and  M  different  units  in  the  basic  alphabet, 
then  at  the  end  of  reaction,  i.e.,  at  equilibrium,  the  total  amount  of  basic  unit) 
will  be  bj .  The  linear  system  Hi  =  b  just  represents  the  fact  that  mass  is 
conserved  in  the  system.  The  matrix  H  is  a  matrix  whose  (i  ,j  )th  entry  gives 
the  number  of  basic  “letter”  j  in  molecular  species  i .  In  the  mathematical 
formulation  it  is  assumed  that  there  are  q  elements  (letters)  in  the  basic 
alphabet  so  that  the  matrix  H  is  q  x  N  and  of  course  the  “state”  vector  (“com¬ 
position”  vector)  is  1  x  N,  where  there  are  np  chemical  species  in  the  pth 
phase  and  N  =n,  +  n2  .  .  .  +  n p.  The  constants  Cyfe  occurring  in  the  expression 
for  F(i)  are  called  “free-energy  constants”  and  are  obtainable  from  tables. 

The  (i,/  )th  component  of  the  state  vector,  namely,  ijj  ,  represents  the  number 
of  moles  of  species  j  in  compartment  (phase)  i .  For  more  detail  on  the  phys¬ 
ical  problem,  applications,  and  extensions  see  Refs  1,  2,  and  7  to  12. 

Now  observe  the  special  form  of  the  matrix  of  second  partial  derivatives 
of  the  function  F(x): 

V2F(x)  =  diagUA^)  4  I,  x^j  I J , 

where  a  is  1,2,.  .  .  ,  0  is  1,  ...  ,p,  and  1(  is  an  N  x  1  vector  with  l’s  in  the 

positions  associated  with  ( i^,  i2j,  .  .  )  as  it  appears  in  the  vector  i  and 

zeros  elsewhere.  The  matrix  of  second  partial  derivatives  of  P  is  thus  given  by 

V2P(x;r)  =  diag(l/*C0  +  2r/xgp)  +  1  h,  [2r/(hjx  -  b,)3]  h,T 


where  hj  is  the  ith  row  of  H,  and  of  course  a  =  1,2,  .  .  .  ,n^;B  =  1,  .  .  .  ,p. 

If  A  =  diag(l/iafl  +  2r/x*/a),  u*  =  h{,  vp  =l[,d(=  2r/(h[i-bj  )3,  and  a\  = 
(  nf  \  p  K 

I  -  1/  E  ijpJ  ,  Eq  2  can  be  put  in  the  form 

V2P(i;r)  =  A  +  .£  ufx\uj  +  Vfttf . 
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Since  A  is  a  diagonal  matrix  and  [a1,],  {cij  3  are  scalars,  it  may  be  seen  im¬ 
mediately  that  Eq  3  is  of  the  desired  form,  and  to  invert  v2P(x;r),  merely  apply 
Eq  1  p  +  q  times.  This  is  not  prohibitive  because  of  the  simple  form  of  the 
matrices  involved. 

The  method  used  in  DeHaven  and  Deland2  as  presented  by  White  et  al  and 
Clasen1’7  is  mentioned  briefly.  This  method  is  in  essence  to  seek  an  approx¬ 
imate  positive  solution  y  =  (yu,yai,  .  .  ,y„  „)T  to  the  system  Hi  =  b.  Then  the 

P * 

free- energy  function  F  is  expanded  in  a  Taylor  series  about  y.  This  can  be 
done  because  y  is  positive.  All  but  the  second- order  terms  are  discarded, 
leaving  a  quadratic  approximation  Q(x)  to  F(x).  This  quadratic  Q(x)  is  then 
minimized  using  a  simple  gradient  descent  method  to  give  a  “better”  approx¬ 
imation  to  the  minimum  of  F.  The  cycle  is  then  repeated  on  Q  using  this  new 
value  for  the  approximate  minimum.  It  is  to  be  noted  that  the  better  the  initial 
estimate  the  more  effective  the  algorithm,  but  convergence  cannot  be  guaran¬ 
teed.  However,  by  applying  SUMT  with  the  knowledge  that  the  problem  is  a 
convex  problem,  convergence  is  guaranteed. 

4.  EXAMPLE 

As  a  numerical  example  of  moderate  size  consider  the  Lung  Level  Res¬ 
piration  Model  presented  by  DeHaven  and  Deland.2  The  data  are  repeated  in 
Fig.  1  with  the  exception  that,  as  indicated  in  Table  1,  an  initial  vector  x  = 

(0,  .  .  .  ,0)  is  used.  In  this  problem,  p  =  7,  q  =  16,  and  N  =  45.  Note  that  in 
DeHaven  and  Deland  11  of  the  56  variables  are  kept  at  zero  and  hence  do  not 
enter  into  the  computations. 

The  solution,  shown  in  Table  1,  was  obtained  on  an  IBM  7044,  32-K, 

IBSYS  version  9.9  computer  in  3.58  min.  This  is  to  be  compared  with  the 
solution  given  by  DeHaven  and  Deland,  where  the  initial  vector  is  somewhat 
“close”  to  the  optimal  vector  and  where  the  method  presented  in  White  et  al1 
was  utilized.  No  computing  times  were  given  by  DeHaven  and  Deland,  and 
hence  a  time  comparison  of  the  methods  cannot  be  made.  It  is  to  be  noted, 
also,  that  there  are  differences  in  the  solution  vectors  for  the  small  com¬ 
ponents.  The  reason  for  this  will  be  invertigated  in  the  future.  Note  also  that 
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TABLE  1 


Comparison  of  DeHaven  and  Deland  Solution  with  SUMT  Solution 


DeHaven  and  Deland2 

SUMT 

Initial  point 

Solution 

Initial  point 

1 

Solution 

6.0000000F-01 

6.4400600E-01 

0.00E-00 

6.4399658E-01 

2.0000000F-01 

2.S895000E-01 

0.00E-00 

2.5896960E-01 

3.0000000E-01 

3.7048000E-00 

0.00E-00 

3.7048090E-00 

3.0000000E-01 

2.9966000F-01 

O.OOE-OO 

2.9966666E-01 

6.7500000F-05 

6.7477199E-05 

0.00E-00 

5.6165829E-05 

7.0000000E-04 

6.9966500E-04 

O.OOE-OO 

6.8800160E-04 

2.1700000E-04 

2.7131800E-04 

O.OOE-OO 

2.061 7855E-04 

2.1000000E-08 

5.3669300E-07 

O.OOE-OO 

1 .1005665E-06 

3.0500000E-07 

8.9100400E-07 

O.OOE-OO 

2.4329989E-06 

5.8500600E-02 

5.7593800E-02 

O.OOE-OO 

5.7145151E-02 

8.0500000E-02 

7.9857600E-02 

O.OOE-OO 

7.9380462E-02 

4.0000000E-03 

3.3371600E-03 

O.OOE-OO 

3.2310024E-03 

2.8950000F.-01 

2.8821 100E-01 

O.OOE-OO 

2.8387379E-01 

1.3500000E-02 

1.4007900E-02 

O.OOE-OO 

1. 38783 36E-02 

9.9200000E-07 

1.4247100E-06 

O.OOE-OO 

3.2831591E-06 

1.8500000E-05 

1.9695800E-05 

O.OOE-OO 

1.7377889E-05 

8.8000000F.-03 

1.1553100E-02 

O.OOE-OO 

1.1551705E-02 

6.8700000F-05 

6.8360400E-0S 

O.OOE-OO 

5.9559742E-05 

4.3400000E-04 

4.2427300E-04 

O.OOE-OO 

4.4194828E-04 

2.2200000F.-04 

2.2238900E-04 

O.OOE-OO 

2.2045284E-04 

2.0900000E-08 

5.5183199E-07 

O.OOE-OO 

1.0947029E-06 

1.1800000E-07 

6.8250600E-07 

O.OOE-OO 

1.8516727E-06 

2.2900000E-02 

2.2458600E-02 

O.OOE-OO 

2.2906 192E-02 

4.6000000F-03 

8.2744299F.-03 

O.OOE-OO 

8.7509019E-03 

5.0000000F-02 

4.4955000E-02 

O.OOE-OO 

4.5060355E-02 

1.7950000F.-01 

1.7888700E-01 

O.OOE-OO 

1 .8322359E-01 

6.1300000F-03 

6.2844200E-03 

O.OOE-OO 

6.3957550E-03 

6.1700000F-07 

1.0429000E-06 

O.OOE-OO 

2.8552963E-06 

5.2800000F-06 

5.6044400E-06 

O.OOE-OO 

7.8058135E-06 

1 .8700000F-02 

2.1 130000E-02 

O.OOE-OO 

2.1 129053E-02 

1 .2000000F-08 

7.973 1300E-06 

O.OOE-OO 

7.4288879E-06 

1.0000000F-07 

3.3106700E-05 

O.OOE-OO 

3.0172473E-05 

5. 8000000 F. -06 

5.5210300E-05 

O.OOE-OO 

5.0560606E-05 

6.4000000F-05 

5.5219900E-05 

O.OOE-OO 

4.8710919E-05 

2.2300000F-03 

2. 1 269900E-03 

O.OOE-OO 

2.141%71E-03 

4.7000000F-05 

1.5P00400E-06 

O.OOE-OO 

2.3368505E-06 

2.6000000F.-03 

1.9890400F.-04 

O.OOE-OO 

1.8211882E-04 

3.6000 OOOE -03 

9.651540OE-05 

O.OOE-OO 

8.5833585E-05 

2.0000 000E-03 

2.8187800E-05 

O.OOE-OO 

2.3550387E-05 

4.5000000E-03 

1  2435600E-03 

O.OOE-OO 

1.2505876E-03 

2.5000000F-03 

7.5461 100E-03 

O.OOE-OO 

7.5730099E-03 

1.0000000E-02 

3.4066600E-04 

O.OOE-OO 

3.0376393E-04 

1.0000000E-02 

4.5441 200E-05 

O.OOE-OO 

3.9017552E-05 

7.7000000E-01 

2.8689100E-02 

O.OOE-OO 

2.8793103E-02 

1.0000000E-01 

I.4953700E-03 

O.OOE-OO 

1.4986067E-03 
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Fig.  1 — Data  for  Lung  Level  Respiration  Model2 

A~  refers  to  an  unoxygenated  1/4  hemoglobin  molecule. 

B"  refers  to  an  oxygenated  1/4  hemoglobin  molecule. 

C~  refers  to  an  amphanion  equivalent  with  -aspect  to  CO2. 

D~  refers  to  an  oxyheme  amphanion  equivalent  with  respect  to  CO2. 
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the  values  of  the  free- energy  function  at  those  minima  yield  essentially  the 
same  minimum.  Neither  solution  satisfies  the  constraints  exactly.  The  point 
to  be  made  is  that  SUMT  requires  no  delicate  analysis  to  obtain  an  initial 
starting  vector.  Also  note  again  that  SUMT  in  conjunction  with  the  rank 
annihilation  of  matrix  inversion  is  extremely  efficient  in  solving  moderate- 
sized  nonlinear  programming  problems. 
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